Scattering a pulse from a chaotic cavity: Transitioning from 
algebraic to exponential decay 

O 

^ ! James A. Hart, Thomas M. Antonsen Jr., Edward Ott, and Steven M. Anlage 

O 

Q 
o 

(N . College Park, MD 20740 



Institute for Research in Electronics and Applied Physics 
University of Maryland 



^ , (Dated: December 20, 2008) 

O 

^ '. Abstract 

The ensemble averaged power scattered in and out of lossless chaotic cavities decays as a power 
' law in time for large times. In the case of a pulse with a finite duration, the power scattered from 

-4— < ' 

' a single realization of a cavity closely tracks the power law ensemble decay initially, but eventually 

S ■ 

I , transitions to an exponential decay. In this paper, we explore the nature of this transition in 

^ ' the case of coupling to a single port. We find that for a given pulse shape, the properties of the 

O ' 

I transition are universal if time is properly normalized. We define the crossover time to be the time 

. at which the deviations from the mean of the refiected power in individual realizations become 

comparable to the mean refiected power. We demonstrate numerically that, for randomly chosen 
■ cavity realizations and given pulse shapes, the probability distribution function of refiected power 

depends only on time, normalized to this crossover time. 
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I. INTRODUCTION 



Waves and wave behavior are ubiquitous. Examples are acoustic waves in matter, elec- 
tromagnetic waves, and physical particles in the quantum mechanical regime. Thus under- 
standing wave behavior is important in many different fields; systems which are radically 
different physically can often be represented by the same mathematics. The simplest model 
of wave behavior is the Helmholtz equation, 

{v' + e)^ = o, (1) 

which typically must be supplemented with boundary conditions. Equation ([T]) describes 
many physical situations exactly (such as acoustic waves within a homogeneous, linear, 
bulk medium or quantum particles in free space). Inhomogeneous situations con often be 
modelled by Eq. ([1]) with k ^ k{r) where k{r) is a function of position. If the system has 
loss or gain, k can be allowed to become complex. Driving terms can be added to represent 
transducers or ports. In this paper, we focus on scalar waves described by Eq. ([T]) with 
constant k, but the results generalize well to many other wave problems. 

Unfortunately, for all but the simplest of geometries, Eq. ([1]) is analytically intractable. 
Thus techniques, both numerical and theoretical, have been developed to solve Eq. ([1]). 
These techniques and their effectiveness vary depending on the regime and physical scenario 
one wishes to study. In this paper, we limit ourself to the semiclassical regime; i.e., the 
regime in which the wavelength of the waves excited in the system is much shorter than 
the scattering elements in the system. In this limit, it is known (via the correspondence 
principle from Quantum Mechanics) that the resulting dynamics are closely related to the 
trajectories a classical particle would take through the system. This analogy applies even 
to purely classical waves, such as waves on the surface of water where the role of classical 
particle dynamics is now replaced by the dynamical evolution of ray trajectories. In this 
paper, we consider only those systems in which the corresponding classical dynamics is 
purely chaotic (i.e., all classical trajectories which start infinitesimally far apart diverge 
exponentially in time). In addition, we focus on the scattering properties of such systems, 
assuming that the system of interest is a closed cavity that couples to the outside world only 
via well-defined localized channels. 



The scattering properties of such wave systems have been well studied, both experimen- 
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Gj and theoretically [7|, Isl, l9|. 



Id, 



11 



, [12, Q, in a wide variety of contexts. 



Much of the the theory has focused on the frequency domain, and sophisticated techni ques 
exist to analyze and characterize the scattering process. See Refs. Is] and 0, Q, Q, Q, Q 



and the references cited therein. Similarly, the time domain response of typical wave systems 
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13| . especially in relation- 



to a delta-function impulse has also been considered 1^, 
ship to fidelity decay(for an overview of fidelity decay, see the Ref. l4] and the references 
therein). In this paper we consider an intermediate situation: we excite the wave system, 
through an external port with a pulse modulated sinusoidal signal, exciting a large but finite 
number of modes. The problem of scattering pulse-modulated sinusoidal waves arises in a 
host of diagnostic situations, such as radar, sonar, nuclear scattering, etc. In what follows, 
for specificity, we discuss our problem in the context of electromagnetic waves. For simplic- 
ity, we consider only lossless two-dimensional microwave cavities excited through a small 
antenna. We emphasize that the results we obtain can be generalized to higher-dimensional 
systems and to quantum mechanical or other wave-chaotic systems(e.g., acoustic or elastic 
wave systems). 

On a formal level, the time domain dynamics of such a system is straightforward. The 
system is open and linear. An incident pulse with a small but finite width in the time 
domain excites a large number of modes in the cavity, which then radiate their energy back 
out through the port. Because the system is linear, the reflected voltage can be expressed 
as a superposition of contributions from modes of the open system. The chaotic dynamics is 
expressed, not through the dynamics of the individual modes, but rather in the eigenvalue 
statistics [l^ and the statistics of the coupling between the port and the cavity. 

As showed in Sec. |Tll the contribution from each mode decays exponentially in time. 
For short times compared with the Heisenberg time (the inverse of the mean spacing of 
mode frequencies), the resulting dynamics will be determined primarily by the semiclassical 
dy,>a,„,cs w,th,„ the cav,ty Q. However, for large times compared with the Heisenben, 
time, the ensemble average of the reflected power decreases as a power law in time |10| . 
This is due to the fact that there is a probability distribution of mode decay rates which 
extends to zero decay rate, and for long times the average is dominated by modes with 
very small decay rates. In the case of a single realization of the chaotic cavity, the incident 
pulse excites a large number of modes with very similar amplitudes, and consequently the 
reflected power initially behaves as though the sum of modes were an ensemble average, and 
the total power decays as a power law. We call this behavior self-averaging. In a single 



specific realization, however, there are only a finite number of modes excited. Eventually 
the slowest-decaying mode in the realization will be much larger than the other modes, and 
the sum will be dominated by this slowest mode, which decays exponentially. Thus for 
extremely long times we expect that the refiected power for any single realization will fall 
exponentially, eventually becoming much smaller than the ensemble average. 

To test this hypothesis, we have created a program that models the time-domain behavior 
of generic chaotic systems. It does this by first generating the spectrum and coupling 
constants of a cavity using the previously published [l7| Random Coupling Model (RCM) 
and then integrating the evolution equations for fields in the cavity, which are modelled in 
the RCM as a set of driven, damped coupled harmonic oscillators. Single realizations of the 
power refiected from these cavities, as well as the ensemble average of 50 different cavities, 
are shown in Fig. [H where we show two very different realizations: one (Fig. [11(a)) in which 
the self- averaging persists throughout the length of the time shown and one (Fig. [11(b)) in 
which self- averaging occurs early, but becomes dominated by solitary slowly decaying modes 
before the conclusion of the numerical simulation. 

Our goal in this paper is to quantitatively describe the transition from self-averaging 
to exponential decay. In particular, we wish to predict the time-scale needed to see this 
transition. In Sec. [Til describe the time-domain model we use for our analysis. In 
Sec. Illlt we find the probability distribution function of the decay rates of the open-cavity 
modes (for the slowest decaying modes in the cavity) as a function of the cavity's port 
refiection coefficient. In Sec. llVl we find the average, standard deviation and (indirectly) the 
higher-order moments of the refiected power as a function of time, and use these moments 
to derive a normalized time which, along with the power spectrum of the incident pulse, 
is all that is needed to obtain a characterization of the transition from self-averaging to 
exponential decay. In Sec. [Vl we evaluate the theory from Sec. [IVl by numerically finding 
the number of modes which fall below certain fractions of the average, and we compare the 
theory with simulation results. 



II. MODEL 



We base our model system on that used in previous work [l7|; specifically a quasi-two- 
dimensional, electromagnetic cavity defined by two conducting plates of area A separated 
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FIG. 1: Using the Random Coupling Model (RCM), we created a program capable of simulating 
the time-domain response of an individual chaotic cavity to a pulse injected into the cavity through 
a small antenna. By repeatedly creating individual cavities using the RCM, we created an ensemble 
of such cavities. The gray lines represent the power reflected back into the cavity from two single 
realizations of the chaotic cavity. The dark solid line represents the reflected power averaged over 
50 realizations of the chaotic cavity. The dashed line represents the time-averaged power for the 
single realization. Figure (a) represents a cavity where self-averaging persists throughout the entire 
simulation, but figure (b) is dominated by solitary modes after about 10~^ seconds. 

by a distance h which are electrically connected along their perimeters by a conducting 
side- wall. The cavity is excited by an antenna that induces currents in the plates. The wave 
equation for this system is 

^^V.,_VX = /.,4 (2) 

where c = (e/i)~^/^ is the speed of propagation of waves in the uniform medium inside 
the cavity, e and /i are the permittivity and permeability of this (non-dispersive) medium, 
VT{x,y) is the voltage difference between the plates, an antenna is modelled through the 
function u{x,y) which gives the profile of current fiowing in the antenna between the surfaces 
(J J dx dy u{x,y) = 1), and I{t) is the time-dependent current driving the antenna. Further, 
as the side walls of the cavity are conducting, Vr — along the perimeter of the cavity. A 
voltage V{t) is induced at the terminals of the model antenna which is given in terms of the 



antenna profile u and Vr 

V = j dxdyuVT. (3) 

The antenna is excited by an incident voltage pulse \4nc(^) arriving along a transmission 
line of characteristic impedance Zq. The incident wave excites the cavity and produces a 
refiected wave pulse Vrc{{t) travelling away from the cavity in the transmission line. At 
the junction between the transmission line and the cavity the voltages and currents at the 
antenna and on the transmission line match, 

V{t) = Mnc(t) + Kef(t), (4) 

I{t) = Zo'[VUt)~Vref{t)]. (5) 

We now introduce Fourier transforms with transform frequency u such that each time- 
dependent variable is represented in the following way, 

VTix,y,t) = I ^e^-^VTix,y,co). (6) 

The transformed field within the cavity is then represented as a superposition of the or- 
thonormal modes of the closed cavity, 

Vt{x, y,Uj) = ^ Cn{uj)(pn{x, y). (7) 
n 

where (V^ ,^ + A;^)0n = 0, and 0„, = on the cavity side walls. 

Solving the transformed wave equation gives the amplitudes Cn{oj) which can then be 
inserted in Eq. ([3]) to find the transformed voltage, 

V{uj) = I{uj)Z,{uj), (8) 

where 



Z,{lo) = -j^J2-^^ Jdxdyu(j)n 



2 

(9) 



is the (exact) cavity impedance. Here k"^ are the eigenvalues of the closed cavity and k = u/c. 



In Ref. [17|, Eq. 14], it was shown that, if one assumed for the purpose of evaluating Eq. iQ 
that the eigenfunctions behave as if they were a superposition of random plane waves, the 
overlap between the eigenfunctions and antenna current profile could be expressed in terms 



of the radiation resistance of the antenna, 



kh fjb f dO ^_^j'\\2 



Rra,{k) = -J^^ / g\m\\ (10) 



where u{k) is the spatial Fourier transform of the profile function u{x,y), and the integral 
is over the angle 6 of the vector k. 

Here Rnad = Re[Zijad] where Zjiad, the radiation impedance, is the impedance V{uj)/I{u) 
that would apply if the cavity side walls were moved to infinity and outward propagating 
radiation conditions were imposed. 

With this random plane wave assumption, the exact impedance Z^, in Eq. was replaced 
by a statistical model impedance, 

7r \ J k^wl RRad{kn) f^^. 

where tf„ are zero mean, unit variance, independent Gaussian random variables. It was 



further assumed in Ref. 



17l | that the eigenvalues fc^ have the statistical properties of eigen- 



values of a Gaussian Orthogonal Ensemble (GOE) random matrix with mean spacing given 
by Weyl's formula, 

(^n+i - kl)n = A = A^/A. (12) 

We now use the relationship (Eq. ([H])) between the voltage Viu) and current lioj) along 
with the transformed version of Eqs. (jl]) and ^ to find the transform of the reflected voltage 
pulse, 

Vreiil^) = p(u;)V^nc(^), (13) 

where the reflection coefficient p{uj) is given by 

Although the derivation above has focused on the electromagnetic case, the expression 
Eq. ( fT4l) describes the reflection of a wide variety of waves when they hit an interface, 
viz., electromagnetic, acoustic, quantum mechanical, etc. The connection becomes closer 
when one considers, as we will, incident pulses whose transformed bandwidth cj^ is narrow 
enough that the radiation resistance and mean frequency spacing can be considered constant 
over the range of excited frequencies. 

The time- dependence of the reflected pulse can be found by using the inverse Fourier 
transformation, 

V^efit) = I ^p(^)l^nc(^)e^'"*. (15) 



The long-term behavior of the reflected pulse is governed by the poles of p{u!) (denoted Uk), 
which satisfy 

Zo + ZM = 0. (16) 

The complex frequencies Uk have positive imaginary parts as they correspond to decaying 
modes. We can approximate the long time dependence of the reflected pulse by pushing the 
inversion contour in Eq. (ITSl) up into the upper half of the cu-plane and deforming it around 
each pole 

where Z'{oj]^ = dZ/duj\^=^i^. Thus, the long time behavior of Kef(^) is determined by the 
properties of eigenfrequencies Uk of the open system. These eigenfrequencies have real values 
whose average spacing is denoted by Auj. In principle, Auj can vary as a function of mode 
number. If we assume that the incident pulse has a spectrum centered at a carrier frequency 
ujq, with a bandwidth tu^ ^ cuq we can relate Au to the mean spacing A of k"^ values 

c^A 

Alu = (18) 

The inverse of this quantity can be identified with what is known as the Heisenberg time in 
the Quantum Chaos community. 

Each mode has a decay rate 7^ = Im(co'fc) which varies from mode to mode. We denote 
the probability density function of these decay rates by ^7(7). Considering the number 
of excited modes to be effectively finite, since each mode decays exponentially, the long 
time behavior of the reflected signal is dominated by modes with the smallest values of 
7fc. From Eq. (fT6!) . along with the expression for Z(uj) in Eq. (fTTl) . it can be seen that 
these weakly coupled modes will have particularly small Wn and thus Heiuk) — knC. Given 
this observation, we can approximate the complex mode frequencies uo^ by solving for the 
poles in the weak coupling approximation. Specifically, in Eq. (fTT]) . our expression for the 
impedance, we separate the term with ujn — knC from the others, 

y( \ v . RRadMAujwl 

Z{uJn)=jXn-J ^, (19) 

where we have changed our indexing labels from kton (because every kn has a corresponding 
ujn), and 

IT ^ k^- kl, kn' 



Thus, we can solve Eq. f ll6p approximately for the complex mode frequencies, 
From this we obtain an expression for the decay rate, 

, 2 RRadZo , . 

^" = '"""'" .(zg + xg - '''' 

The reactance X„, like the impedance Z is a statistical quantity. It has an average value 
to which all the terms in Eq. (1201) contribute, and which can be calculated by replacing the 
sum by an integral |17i |. 

iX„,^X.^^-lp{fK^^^}. (23) 

where the symbol P indicates that principal value definition of the the integral is to be 
taken. This average value is the radiation reactance of the antenna. The reactance X„ has 
a fluctuating part which scales as the radiation resistance and is due primarily to terms in 
the sum where n and n' are not too different, 

Xn = Xjiad + RRadin- (24) 

The quantity ^„ has a universal distribution which we will investigate in depth later. 

Using Eqs. (1191) and (l^Tl) we may evaluate Z'{uJn) in the denominator of Eq. (IT7|) . The 
result for the reflected signal is 

Kef(t) = -25^ ^f^^' wle^^-'AuVU^^). (25) 
Taking the magnitude of this, we obtain the reflected power, 

Pref(t) =Pref(t)+Pref(t), (26) 

where 

I ^ 1 2 

PUt) = J2 ^^±Jl!L^e-"^-\ (27a) 

and ipn is the phase of Zq + jXn. 



The two contributions to the reflected power fl27ap and fl27bl) are very different. In the 
first contribution the terms decay exponentially and smoothly and the sum is always positive. 
In fact, if we smooth over a timescale longer than the Heisenberg time, this first term will 
remain essentially unchanged. The second term, on the other hand, oscillates rapidly on 
a timescale comparable to the Heisenberg time, but tends to zero if averaged over long 
timescales. For the very long timescales needed to see the transition from self-averaging to 
exponential decay, we can treat the rapidly fluctuating terms in Pref(^) as random variables 
with the phases in the exponents ((cu^ — cjj^) t) being uniformly distributed. Under this 
assumption, we find that, for a single realization of the chaotic cavity, the fluctuating part 
of Pref is random and has a variance of 

a2=([P,ef(t)]\<PrefW. (28) 

where {■ ■ ■)t indicates a sliding averaging in t over a timescale that is long compared to the 
Heisenberg time but short compared to the characteristic time for variation of Prcf (t) . That 
is, the order of magnitude of the oscillating part of Prof is typically the same as that of 
the smoothed part. Thus, if the smoothed part of Prcf drops exponentially, the fluctuations 
around it will as well. Hence, if the power stays self-averaged, the fluctuations will be as 
large as the signal itself. When we consider the transition from self-averaging to exponential 
decay, we consider only the statistics of the smoothed part of P^ef, ignoring the oscillating 
part which does not contribute to the self-averaging. Thus in our theory we consider only 
the time-averaged power Pref(t), Eq. f l27ap . which is the key result of this section. 



III. FINDING Pji'jn) 

From Eq. (12 Taj) , we see that the average reflected power is a sum over contributions from 
exponentially decaying modes. Because of the exponential decay, the relative amplitudes 
of the modes will separate exponentially in time, with the modes with the smallest 7„ 
eventually dominating the sum. Thus, the crossover time from self-averaging to exponential 
decay depends on the behavior of the probability distribution function of 7„ for small values 
of 7„. In this section we find the behavior of P^(7„), the probability distribution function 



for the decay rates for 7„ ^ Au. Previous work 



nas 



in the case of a lasing chaotic cavity, see Refs. 18 



3een done on the subject (for instance. 



19|). including analytical solutions for 



the P-yi'y) for all 7 20|, |2l|], but because we focus on the single port case with time reversal 
symmetry for small 7 only, many approximations can be made which greatly simplify the 
derivation, which we present here. 

We start by considering the statistics of C,n, where is defined in Eq. 1^^. We show in 
Appendix El that the statistics of ^„ are given in terms of the angle ipn = tan~^(^„), where 
ipn is distributed according to the pdf, 

= (29) 

Using this result and Eq. (!22|) . we find an expression for P^(7n) where 7^ -C Aa;: 



Pliln) = / #n COsV'n / dwC ^'^Si 7„ 



w 



27r 7-^/2 Jo \ vr [1 + (r^ tan(?/'„) + Xr^nY] 



■n/2 

(30) 

where = Ruadik) / Zq and = X^adik) / Zq. The innermost integral can be evaluated 
leaving only an integral over ipn- Further, since we are only interested in the case of small 
7„ ^ Ac(j, the main contribution comes from \w\ <^ 1. The result is 

P,(7„) = —2= for 7„ « Au;, (31) 
2V7nAa; 

where 

Po = (2r,.)"^/^ / 1/ cos^ + {rr svciipn + cos?/'„)2. (32) 

The quantity Pq given in Eq. ( l32l) can be rewritten in terms of the radiation reflection 
coefficient of the port that applies when the walls of the cavity have been moved out to 
infinity, 

= Iti' ^^^^ 

where Zr = Tr + ixr = {Rnad + jXRad) / Zq is the normalized radiation impedance of the 
antenna. To see this, we introduce the intermediate variable (3 = — 1 and define a new 
integration variable (p = ipn — arg(/3)/2 in Eq. ( l32i) . The result of these variable changes is 



^o = l/2^^i?(^^|, (34) 



\Pr\ V i - \Pr 



where 

"77/2 



E{k) = J d(f)^l~k^sm\(f)) (35) 
is the complete elliptic integral of the second kind. 
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FIG. 2: A comparison of numerically generated values for Pq (circles) with the theoretical result 
from Eq. (|34p (the solid line). The circles represent numerical calculations of Pq with the radiation 
reactance of the port set to be Xjiad = 0. To get different values of Ip^l, was changed as 
described in Eq. ([36|) . 

We confirm Eqs. fl3ip and fl34p numerically by generating an ensemble of 7„ values. To do 
this we solve Eq. 0161) by generating different realizations of the Gaussian random variables 
Wn and random matrix eigenvalues k"^ appearing in the definition of Z{uj), Eq. (ITTl) . We 
find the mode frequencies by noting that a.s Zq —>■ oo, Un ^ knC for all modes. We then 
introduce Yq = Zq^ and differentiate both sides of Eq. (fT6|) with respect to Yq, obtaining a 
differential equation for (^^(Fo), 

dYo ~ Z'M ' ^^^^ 
which can be solved numerically to find uJn for finite Zq. Note that although both Z^{ujn) 
and Z'{uJn) are singular as Un — > knC, their ratio is finite. 

By generating 1000 different realizations of k'^ and uj^ (truncating the spectrum to include 
only 600 terms), and integrating Eq. (136|) numerically using fourth-order Runga-Kutta from 
lo = to lo = -R/jad' is possible to generate pdfs of Wn = as a function of \pr\- We 
choose the pdfs of Wn instead of 7„ because Pu,{w = 0) = Pq/Au, which is finite and thus 
numerically easier to fit. The results are shown in Fig. [2] where the numerical results and the 
theory are seen to be in clear agreement. We note that this numerical test (solving Eq. (|36|) 
for Yq = R~j^\^ does not assume the weak coupling limit and thus confirms our assumptions 
in obtaining Eq. fl34l) . 



IV. THE STATISTICS OF Pref(i) 



The smoothed reflected power -Pref(t) given by Eq. fl27ap is a sum of terms each of which 
is a random variable. The terms are not strictly independent. This follow from the fact that 
there are correlations between the eigenvalues of the closed system, and 7„, given by Eq. (I22l) . 
depends on these eigenvalues through the reactance X„, deflned in Eq. fl20l) . Fortunately 
the correlation is signiflcant only for almost adjacent modes. For times large enough that 
the self- averaging breaks, the fraction of modes contributing will be small, and thus, the 
majority of contributing modes will be well separated and approximately independent of 
each other. 

Hence for our purposes, Pref can be treated as a sum of a large number of independent 
terms. Thus, for times when a large number (but small fraction) of modes have compara- 
ble magnitudes, for an ensemble of cavity realizations, Pref is a Gaussian random variable 
centered on (Prof(t)) with a small standard deviation. As we demonstrate in the following 
sections, the standard deviation starts out small, but as the number of contributing modes 
decreases, the standard deviation increases relative to the mean, eventually becoming much 
larger than the mean. As this happens, the simple Gaussian distribution changes into a more 
complex distribution with the majority of modes becoming much smaller than the average, 
corresponding to the shift from self-averaging to exponential decay. 

These shifts can be treated analytically by considering the moments of Pref. We first 
(Sec. IIV A|) consider the mean and standard deviation of Pref to find a scaling law de- 
scribing the transition from Gaussian to non-Gaussian behavior. Armed with the results 
from this comparison, in Sec. IIVBI we generalize the results to higher-order moments (via 
the cumulants), showing that for large times all moments of Pref obey the same scaling law. 
We then numerically demonstrate that the cumulative distribution function of Pref/ (Pref) 
satisfies the scaling law for multiple pulse shapes, as predicted. 



A. The Mean and Variance 

We can calculate the mean and the variance of Pref for all times as 

/p \ >^ |27rAa;1^rie(^n)P 

(^ref) = 2^ y /^l, (37) 



and 







(39) 



where 

Evaluation of the integral in Eq. fl5^ gives 

Equations (1371) and ( l39l) give the result that the average reflected power (averaged over an 
ensemble of reflecting cavities) decreases as a power law in time, which is in agreement with 



previous theory 



10 



(P,ef(t))~t-'/'. (41) 



Equation fl38l) is useful for finding the range of values that are most likely to contain Pref; 
for small times with an approximately Gaussian pdf for Pref, we expect that the majority of 
realizations will fall within the range [(Prcf)— 2crp, (Pref)+2crp] where ap = ((Prcf— (Prcf))^)^^^- 
For large times, however, ap > (Prcf)- We see this by first considering the ratio 

^ {Aujty/' r(9/2) 
/if Po 27/2r(5/2)2- ^ 

Thus, for large times, /i2 ^ yU?, and /i2 dominates Eq. fl38l) . For large times, we have 

_al_ ^ {tAuj f' r(9/2) Enl^nc(a^n)r .4^^ 

Equation (143!) can be made more transparent by considering the sums over iV^ncP™- The 
incident pulse can be considered to have two independent properties: a shape and a width. 
If we double the width of the pulse in the frequency domain (or equivalently if we halve the 
average mode separation) without changing the shape, the sums in Eq. ( H3l) will, to a good 
approximation, simply double. We thus define the effective number of modes excited by the 
wave to be 

^ ^ En|V^nc(^n)|f . . 

EJ^nc(c.n)r ■ ^ ^ 

In the case of a square wave excitation in the frequency domain, Eq. (j44|) gives exactly 
the number of modes excited. In the case of more typical pulses, such as a Gaussian pulse, 
Eq. (H4l) defines a relationship between the pulse width and the number of significant excited 
modes. 



Substituting Eq. (jB]) into Eq. (jl3]), we get 



{F„f)2 2'/2r(5/2)2 

where 

r = (46) 

As long as ap/ (Pref) is small, it is reasonable to expect the majority of realizations of Pref to 
be within two sigma of the average, and numerically we find that this is true. From Eq. fHSl) . 
we see that for tAu ^ 1 and r ^ 1 (possible because is assumed to be large) this is 
possible. Eventually the standard deviation will be comparable to the mean and for very 
long times the standard deviation will be much larger than the mean. This shift corresponds 
to the change from self-averaging to exponential decay. 



B. Higher Moments 

An analysis of the higher moments of Pref follows essentially the same steps as those to 
find the mean and variance. We find the moments of P^ef by finding the moments of the 
individual terms in Pref, dropping all but the leading order term in t"^/^, and combining 
them properly to get the moments of the sum. We cannot do this by simply summing the 
moments of the individual terms; the sums of the moments are not in general the moments 
of the sum. However, if we define the moment-generating function, 

M{h) = (e'^^-O = 1 + Y1 , > (47) 

we see that the moments of Pref are given by 

(0 = M(")(0). (48) 

Here M^"^^{h) is the mth derivative of M{h) with respect to its argument. This can be 
related to a function known as the cumulant-generating function 

°° hV 

g{h)=\og{M{h)) = J2^P^ (49) 



where Km is the mth cumulant, defined as 



Km = 9^"'\0). (50) 



We show in Appendix ([B]) that, in analogy to Eq. fl43p . the higher-order cumulants (and 
thus all higher-order moments) of Pref are given by 



r(2m+l/2) AT'-i^J 






Mnc(c^n)P)''" 



/ ^2m+l/2F('c; /OW I'Y^ ITZ (,, \\2\'<^ ' ' 



If we use the definition of from Eq. (jH]) and approximate all sums over n with integrals 
over ojn, we find that the expression N"^~^ J2n l^nc(^^n)P™'/ {J2n l^nc(^n)n™ is, to a good 
approximation, independent of the width of the power spectrum but dependent on the shape. 
In the case of a square power spectrum, this factor is identically one for all m. For a Gaussian 
pulse we find that 

(52) 
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Similarly, for a pulse with a Lorentzian power spectrum. 
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(53) 

Equation ( 15T1) . combined with replacing the sums over |Mnc(^^)P™ with integrals, demon- 
strates the most important theoretical result of this paper: all statistical properties of the 
reflected power depend only on the shape of the pulse (independent of width) and the nor- 
malized time T defined in Eq. fj46l) . Thus the cross-over from self-averaging to exponential 
decay, no matter how measured, will depend only on r and the pulse shape. 



V. NUMERICAL RESULTS 

In this section, we compare different methods of calculating -Pref(i^) to show that our 
theoretical conclusions are correct. To view the resulting distributions, we find the ensemble 
average of the calculated values of Pre{{t) and then compare the individual realizations to 
the average. In particular, we define C{a, r) to be the fraction of realizations which are less 
than a times the ensemble average (i.e. C{a, r) is the cumulative distribution of Pref at the 
normalized time r). 

To both test and evaluate the theoretical results in Sec. IIVI we perform two separate, 
independent calculations which should, according to our theory, produce the same results. 
The first method calculates the sum in Eq. f l27al) with the 7„ independent of Re(ci;n) and 
distributed according to the Porter-Thomas distribution with one degree of freedom, 

p-7/2 



This distribution is chosen because it has the same behavior for small 7 as is indicated in 
Eq. fl3T|) . We consider two different pulse spectra, Mnc(i^n), Gaussian and Lorentzian, with 
two different widths = 20 and 30, where is defined in Eq. dH]). Finally, we take the ujn to 
be uniformly spaced when evaluating the sums. We call these results the theoretical results 
because they are a numerical evaluation of the theoretical assumptions used in Sec. IIVI The 
theoretical results are shown in Fig. [3] for the case of the two pulse shapes and two spectral 
widths. The first thing to note about the plots is that the results for = 20 and = 30 
lie on top of each other, showing that the definition of r (l46!) correctly accounts for variation 
of the pulse width. (There is a small deviation in the Lorentzian case for small values of 
a that will be addressed subsequently.) The second thing to note is that the C{a > 0.3) 
curves for the two pulse shapes are very similar. Thus, the fraction of realizations close to 
or greater than the mean is the same in the two cases. Where the two pulse shapes differ 
is for times r ^ 1 and small a <^ 1. In the Gaussian case almost all realizations fall well 
below the average as r gets large, whereas in the Lorentzian case there is a larger fraction 
of realizations with measurable power (a > 0.001) at late time. This is due to the long tail 
in the Lorentzian distribution exciting a large number of modes with small but significant 
levels of power. The difference between the A^ = 20 and A^ = 30 cases is due to truncation 
of the spectrum at 600 modes. 

The second test employs the time-domain code used to generate the data in Fig. [H 
We then time-smooth the resulting power (using a Gaussian window with a width of 10 
Heisenberg times) to calculate Pfcf- The time domain code is described in Appendix O In 
Fig. (jl]) we compare results for C(a,r) using 50 realizations with the theoretical curves. 
The time-domain code is run only to r = 1 which for these parameters corresponds to 1744 
Heisenberg times. The time domain simulation results agree well with the theoretical results 
considering the finite sample size. 

In addition, we have performed tests which have allowed the value of Pq fo vary, and have 
solved Eq. (136|) to get the complex values of ui^- The results agree well with the theoretical 
results of Fig. [3] and are not displayed. 
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FIG. 3: The fraction of realizations of Pref which are less than a(Pi.ef) ^ a function of normalized 
time r for (a) a Gaussian spectrum and (b) a Lorentzian spectrum. The black lines('+' symbols) 
represent the statistics for N = 20(30). Note that plots for = 20 and = 30 are slightly 
different for the Lorentzian case with small a. This is due to the fact that the contributions for 
small a come from the tails of the distribution, which we numerically truncated to calculate these 
plots. 
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FIG. 4: The fraction of realizations of Pref which are less than a(Pi.ef) ^ a function of normalized 
time r for the theoretical results calculated numerically (the solid lines) and the same results 
calculated from integrating Eq. ([2]) directly (indicated by the '+' symbols). Random Matrix Theory 
is explicitly used to calculate the spectrum and coupling constants for the time-domain integration. 



VI. CONCLUSIONS 



In this paper, we have found numerically and theoretically that the long term behavior of 
power reflected from a lossless, microwave cavity excited through a single port self-averages 
for times larger than the Heisenberg time, decaying as a power law in time. We have also 
found, theoretically and numerically, that for times much longer than the Heisenberg time, 
when r, the normalized time, is of order 1, that single modes in the cavity will begin to 
dominate the long term decay and the reflected power will begin to decay exponentially. The 
details of this behavior have been found to depend on the shape of the power spectrum of the 
incident pulse that excited the cavity, but to otherwise depend only on the normalized time. 
Because much of the theory used to derive this behavior depends only on generic Random 
Matrix Theory, we expect that this behavior will translate into other lossless wave-chaotic 
systems (e.g., acoustic, quantum mechanical, etc.), independent of details. 
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APPENDIX A: FINDING THE DISTRIBUTION OF ^„ FOR SMALL 7„ 

To find the distribution of defined in Eq. ( I24l) for small 7„, we exploit the fact that, 
in a two-port system with the ports identical and described by Random Matrix Theory, the 
diagonal elements of the normalized impedance matrix each have the same statistics as the 
single-port normalized impedance. Then using the exact statistics of the two-port RMT 
impedance, we can find the statistics of the one-port impedance ( l20l) . 

We see this by first writing the elements of the two-port normalized impedance matrix 
as a sum, analogous to Eq. (ITT]) . 



where the independent Gaussian random variables and the k"^ have the statistics of 

the eigenvalues of a GOE random matrix. 




(Al) 



sm 



(A2) 



As shown in previous work [23], the 2x2 matrix ^ has the following statistics: its eigen- 
values tan6'i, and tan ^2 have a joint pdf, 

2 

and its eigenvectors (cos z/, sin z/) and (— sin z/, cos z/) have v uniformly distributed and in- 
dependent of 9 1 and 62. Consequently, a diagonal element of ^ can also be parameterized 

as 

^i^i = cos^ r] tan 9i + sin^rj tan 62- (A3) 

Comparing Eqs. (]A1|) and (]A3p . we see that the singularity at A; = A;„ in Eq. (lAip 
is matched by either di or 62 going through 7r/2; for specificity we assume that it is 9i. 
For small 7„, corresponding to small w^, the coefficient of the singularity is small, which 
corresponds to cos^r] ^ 0. Thus, for small 7^, in has the statistics given by 

= tan6'2|ei=7r/2 (A4) 
which inserted into Eq. (]A2|) produces the pdf for ipn = tan~^ ^„ = 62 

P(;M = ^ (A5) 

Numerically we confirm this by generating a single 600x600 element matrix from the 
Gaussian Orthogonal Ensemble and calculating and scaling the eigenvalues to get an ap- 
propriate spectrum. We then repeatedly generate 600 realizations of 600 coupling constants 
and use them to calculate 360,000 realizations of X„, which we then normalize to calculate 
ipn- The resulting statistics are demonstrated in Fig. [51 

APPENDIX B: FINDING THE CUMULANTS OF Pref 



To obtain Eq. (15T|) . we note that the cumulant generating function of Prcf obeys 



This result is a specific example of a general property of cumulants 



(Bl) 



24| : The mth cumulant 



of a sum of independent variables is the sum of the mth cumulants of the single variables. 
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FIG. 5: A comparison of a numerically-generated pdf of ipn ('+' symbols) with the anticipated 
result from Eq. (|A5p . cos('(/'„)/2 (the solid line). 

Thus, in analogy to Eq. fj49l) . we define the cumulant-generating function and the cumulants 
kp for each term in the sum in Eq. (IBip as 



7n ^-2'y„t 



g{q) = log ( ( exp ( q^^e 



and by matching coefficients of h"^ in Eq. (IBll) . we get that 



p=l 



(B2) 
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All that is left is to find the long-term behavior for km- To do this, we note that we can 
rewrite the average exponential in Eq. flB2l) as 



1: 



n -2-ynt 



q fir. 
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Because g is a dummy variable which can be arbitrarily small, we can also expand the 
logarithm in Eq. ( ]B2p to get that 



q'^f^n 
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By matching coefficients of on both sides of Eq. (]B5|) . we find that 24 1 



Ki = /ii. 



= /i2 - /ii. 
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-1)>^ 



(B5) 

(B6) 
(B7) 
(B8) 



where the ehded terms are products of different /i„ such that the indices add up to m. For 
large t, all of these polynomial terms are small compared We can see this by noting 
that //^ oc (tAu;) -2^-1/2 _ xhus 

oc (^A^)l/^ (B9) 

where the proportionality constant can be shown to be order 1. For every extra factor of 
Hi included in a term, we pick up an extra factor of (tAcj)^/^ in the numerator of the ratio 
between /i^ and that term. Thus for large times we have that /x^ is much greater than any 
of the other polynomial terms in Eq. (1B8I) and therefore 

K„ ^ /im. (BIO) 

Combining Eqs. (gO]), ([63]), and fiBTOl) . we get Eq. 



APPENDIX C: THE TIME DOMAIN CODE 



In this section, we describe the time domain code used to create the realizations in Fig. [H 
This code effectively solves Eq. ([2]) using the approximations that were inserted into Eq. ([9]) 
to produce Eq. ( fTTl) . In addition, it makes use of a slowly varying envelope approximation 
which greatly increases the size of the numerically stable time-step and also transforms 
Maxwell's Equations into Schrodinger's Equation. 

To solve Eq. ([2]), we first find expand V{x, y, t) in terms of the eigenfunctions of the closed 
system 



Cn{t)(j)n{x,y) 



(CI) 

' J de \uiwoc)\^ 

We note that the Cn{uj) from Eq. ([7]) are proportional to Fourier transforms of the c„(t). 
Substituting Eq. into Eq. and using the orthonormality of the 0„, we get 



1 2~ 

^^^C„(t) + k^Cn{t) 



87iRr{ujoc) dl{t) 



UJQ 



dt 



J dx dyu(j)n 
f dO \u{wqc)\' 



(C2) 



where we have used the definition of radiation resistance from Ref. 17|, Eq. 19] to remove the 
factor h^. The value of ujq is the modulation frequency used in the envelope approximation 
(See Eqs. and (ICiD l 

To apply the envelope approximation, we assume that 



I{t) = /env(t)e- 



(C3a) 



5n(t) = d^{t)e^^°' 



(C3b) 



where 



, , -^env 



d 

dt 



d^{t) < UJodm{t). 



d"^ d 

-^dmit) < UjQ — dmit). 



(C4a) 
(C4b) 
(C4c) 



Then we drop all terms which are small, noting that A;„ ~ t^o/c, which implies that fc^c^ 
= {knC — ujQ){knC + ojq) is on the order of ujq. This gives us 



2juJod_ 
c2 9t 



dn{t) = 8jTlRR{uJoC)Ienvit) 



J dx dy 



uc 



J dO \u{w{)c) 1^ 



(C5) 



Again we replace the overlap integral between 0„ and u with the statistical approximation 
found in Ref. 



17|, Eq. 14] to get 

+ {K - „2 



dn{t) = V8AWnjRR{uJoC)Ienv{t). 



(C6) 



dt ' " 

Similarly, combining Eqs. (ICip and ([3]) and using the envelope approximation throughout, 

we get 

V;nv(t) = $^K(t), (C7) 

n 

where Knv(^) is the envelope of V(t) in analogy to Eq. fICSP and 



Vn{t) 



4tt 



dn{t)Wr. 



(C8) 



Solving Eqs. (jlj) and ([5]) for I{t) by eliminating Kef(^) and inserting the result into 
Eq. (IC6]), we get 



c2 9t 



where Vi,env(^) is the envelope of Vinc(i) in direct analogy to Eq. ( 1C3I) . 

Equation ( ICOP is a set of complex first order linear differential equations analogous to 
Shrodinger's equation. By truncating the spectrum to a finite number of modes, it is possible 
to solve Eq. (\C9\i numerically via standard numerical integration techniques. In our case, 
we choose forth-order Runga Kutta. We generate the values of /c^ — by generating 



600x600 random matrices from the Gaussian Orthogonal Ensemble, finding the spectrum, 
and unfolding it such that the — have a uniform density. We also generate the 600 
Wn as Gaussian random variables with mean and width 1. All of the remaining variables 
(including the initial conditions) are physical parameters that must be set to match the 
situation we wish to simulate. 

For the runs displayed in this paper, we chose Rr{u!oc)/Zq — 1, cuq — 22.5 GHz, and 
A = 10 m~^. The kn were chosen to lie between Ri 51 m~^ and 93 m~^. For initial conditions, 
Vn{0) — 0. The envelope of the incident pulse, Vi^env, had the form 

V.,Ut) = e-^"^--'^'/' (CIO) 

with (7^ = 150 MHz. 
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